clear all
clc
close all

w = 2*pi / (12*60*60); % Угловая скорость вращения радиус-вектора НС
Rz = 6371e3;
Rsv = 19100e3 + Rz; % Радиус орбиты НС
c = 3e8;
f0 = 1602e6;

Tmod = 6*60*60;
dT = 5*60;
T = 0:dT:Tmod;

rr = [0;Rz];

figure(1)
subplot(2,2,1)
p = plot(0,0,'*',0,Rz,'*');
axis([-Rsv Rsv 0 2*Rsv]);
xlabel('X');
ylabel('Y');


subplot(2,2,3)
p2 = plot(0,0);
axis([0 Tmod 1.5e7 2.8e7]);
xlabel('time');
ylabel('distance');

subplot(2,2,2)
p3 = plot(0,0);
axis([0 Tmod -1e3 1e3]);
grid on;
xlabel('time');
ylabel('Vr');

subplot(2,2,4)
p4 = plot(0,0);
axis([0 Tmod -6e3 6e3]);
grid on;
xlabel('time');
ylabel('fd');

D = nan(1,length(T));
V = nan(1,length(T));

for i = 1:length(T)
    rsv = Rsv * [cos(w*T(i)); sin(w*T(i))];
    set(p(1),'XData', rsv(1),'YData',rsv(2));
    set(p(1),'color','r');
    pause(0.05);
    
    rv = rsv - rr;
    
    D(i) = norm(rv);
    set(p2(1),'XData',T,'YData',D,'LineWidth',2);
    
    if (rsv(2) > Rz)
    V(i) =  (D(i) - D(i-1))/dT;  
    end
    
    set(p3(1),'XData',T,'YData',V,'LineWidth',2);
    set(p4(1),'XData',T,'YData',V*f0/c,'LineWidth',2,'Color','r');
    
end

figure(2)
hist(V*f0/c)

rsv = Rsv * [cos(w*T); sin(w*T)];
rv = rsv - repmat(rr,1,length(T));
x_rv = rv(1,:);
y_rv = rv(2,:);

D = sqrt(x_rv.^2 + y_rv.^2);

V = diff(D)/dT;
fd = V*f0/c;

figure(3)
plot(T(2:end),fd);
